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I.     INTRODUCTION 

Finite  impulse  response  (FIR)  filters  have  seen  widespread  use  in  the  field  of 
two-dimensional  digital  signal  processing.  Reasons  for  this  include  the  inherent 
stability  of  FIR  filters,  the  ability  to  achieve  a  linear  (or  zero)  phase  characteristic 
in  the  frequency  response,  and  relative  ease  of  design.  Traditional  approaches  to 
the  design  of  such  filters  include  the  window  method,  in  which  a  smoothing  window 
is  applied  to  the  Fourier  series  coefficients  of  the  ideal  frequency  response;  the  fre- 
quency transformation  method,  in  which  a  one-dimensional  (1-D)  prototype  filter 
response  is  mapped  into  a  function  of  two  frequency  variables  using  an  appropri- 
ate transformation  function;  and  the  frequency-sampling  method,  in  which  desired 
frequency  response  values  are  specified  at  certain  sample  points  in  the  frequency 
domain.  Of  the  three  methods  described  above,  the  frequency-sampling  method  is 
probably  the  least  utilized  and  the  least  understood,  and  hence  warrants  further 
investigation.  This  thesis  deals  with  a  comparative  analysis  of  various  methods  of 
designing  two-dimensional  (2-D)  FIR  filters  using  frequency-sampling  techniques. 
In  the  analysis,  some  new  methods  for  frequency-sampling  design  are  presented. 

The  2-D  frequency-sampling  design  method  is,  in  essence,  an  interpolation 
problem  in  two  variables  where  the  two  variables  are  the  two  frequencies  of  interest, 
denoted  as  ujx  and  u2.  The  frequency  response  is  specified  at  selected  sample 
points  in  the  {uj1,<jj2)  plane,  and  is  expressed  in  terms  of  linearly  independent 
basis  functions.  The  solution  to  this  interpolation  problem  involves  finding  the 
coefficients  associated  with  the  respective  basis  functions  in  order  to  meet  the 
interpolating  conditions.  These  coefficients  are,  in  general,  closely  related  to  the 
impulse  response.     The  various  approaches  to  the  frequency-sampling  design 


problem  involve  differing  constraints  placed  on  sample  location  in  the  (u^  , oj2)  plane 
and  differing  sets  of  basis  functions.  Because  the  resulting  frequency  response  takes 
on  exactly  the  prescribed  value  at  each  frequency  sample  location,  the  response  is 
somewhat  predictable  at  the  onset  of  the  design  process.  A  disadvantage  of  2- 
D  frequency-sampling  design,  however,  is  that  the  existence  and  uniqueness  of  a 
solution  is  dependent  upon  the  location  of  selected  frequency  samples  and  the  basis 
functions  specified;  in  other  words,  degenerate  cases  may  arise  in  which  no  solution 
exists  to  the  bivariate  interpolation  problem. 

FIR  filter  design  techniques  denoted  as  "frequency-sampling"  traditionally  in- 
volve sampling  a  desired  frequency  response  at  equally-spaced  intervals  and  then 
solving  for  the  impulse  response  coefficients  using  an  inverse  discrete  Fourier  trans- 
form (IDFT)  [Ref.  l].  For  the  2-D  case,  the  process  is  similar,  with  the  desired 
frequency  response  sampled  at  equally  spaced  points  on  a  Cartesian  grid  in  the 
(u!1,u>2)  plane  [Refs.  2,  3].  This  approach  will  be  referred  to  as  the  "uniform 
sampling"  method  in  this  thesis.  There  are  two  significant  benefits  to  such  an  ap- 
proach in  two  dimensions.  First,  with  frequency  samples  placed  at  the  vertices  of 
an  N  x  N  Cartesian  grid,  the  existence  of  a  unique  filter  satisfying  the  interpolating 
conditions  and  possessing  an  impulse  response  with  an  N  x  N  region  of  support  is 
guaranteed.  Second,  use  of  the  IDFT,  or  alternatively,  the  fast  Fourier  transform 
(FFT)  algorithm,  assures  a  computationally  efficient  design  process.  The  princi- 
pal disadvantage  of  the  approach  is  that  there  is  no  flexibility  in  the  placement  of 
frequency  samples.  The  approach  is  presented  in  more  detail  in  Chapter  II. 

A  different  approach  to  frequency-sampling  design  of  2-D  FIR  filters  involves 
allowing  frequency  samples  to  occur  anywhere  in  the  (wj ,  uj2  )  plane  and  fixing 
the  region  of  support  for  the  filter  impulse  response.  This  approach  is  presented 
in  Chapter  III  and  will  be  denoted  as  the  "arbitrary  sampling"  technique.    The 
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method  is  quite  flexible  since  samples  can  be  arranged  to  enhance  certain  desired 
features  in  the  resulting  frequency  response,  such  as  flat  passband,  sharp  roll-off, 
equi-ripple  behavior,  etc.  The  design,  however,  involves  the  solution  of  a  large 
system  of  linear  equations  for  a  filter  of  any  practical  order  and  hence  it  is  quite 
computationally  intensive.  Additionally,  with  no  constraints  on  the  location  of 
samples,  degenerate  cases  may  arise  in  which  no  solution  exists  which  meets  all 
interpolating  conditions  for  the  specified  impulse  response  region  of  support. 

A  novel  2-D  frequency-sampling  design  technique  is  presented  in  Chapter  IV, 
and  is  referred  to  as  the  "nonuniform  sampling"  method.  This  method  places 
certain  constraints  on  frequency  sample  location  but  still  possesses  much  greater 
flexibility  than  the  uniform  sampling  approach.  Rather  than  solving  one  large 
system  of  equations,  this  technique  involves  the  solution  of  several  smaller  systems 
of  linear  equations.  Given  the  constraints  on  the  sampling  geometry,  existence  and 
uniqueness  of  a  solution  is  guaranteed.  Additionally,  the  design  method  is  less 
computationally  intensive  than  the  arbitrary  sampling  case  but  more  intensive,  in 
general,  than  the  uniform  sampling  case. 

In  comparing  and  contrasting  these  various  approaches  to  2-D  frequency- 
sampling  design,  several  interesting  points  arise.  First,  frequency-sampling  design 
techniques  involve  some  sort  of  decision-making  process  in  determining  where  to 
place  the  samples  in  order  to  best  match  the  desired  frequency  response.  Methods 
which  are  more  constrained,  such  as  the  uniform  sampling  approach,  involve  a  lesser 
degree  of  decision-making  and  hence  can  be  considered  easier  to  use  from  that 
perspective.  Second,  by  reducing  the  number  of  constraints  on  the  placement  of 
frequency  samples,  flexibility  in  matching  the  desired  frequency  response  is  gained. 
Finally,  and  perhaps  most  significantly,  computational  efficiency  in  determining 
filter  parameters  can  be  directly  related  to  the  number  of  constraints  on  sample 


location.  Of  the  techniques  discussed  here,  those  which  are  more  constrained  tend 
to  be  more  computationally  efficient.  Thus  there  are  tradeoffs  to  be  considered 
when  selecting  one  of  these  techniques.  In  this  thesis,  a  detailed  analysis  of  these 
concepts  is  performed  and  the  quality  of  these  frequency-sampling  methods  in 
matching  desired  filter  characteristics  is  investigated  through  design  examples. 
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II.     DESIGNS  WITH  UNIFORM  SAMPLING  ON  A  CARTESIAN  GRID 

The  approach  to  frequency-sampling  design  of  2-D  FIR  filters  discussed  in  this 
chapter  is  the  one  referred  to  as  "frequency-sampling"  in  various  texts  [Refs.  1, 
2].  It  will  be  referred  to  the  "uniform  sampling"  method  through  the  remainder  of 
this  thesis.  This  method  has  significant  computational  advantages  through  use  of 
an  inverse  discrete  Fourier  transform  (IDFT)  approach.  Additionally,  this  method 
will  never  result  in  any  degenerate  cases,  i.e.  situations  in  which  a  filter  satisfying 
each  of  the  interpolating  conditions  cannot  be  found.  However,  due  to  constraints 
on  the  location  of  samples,  the  method  is  rather  inflexible. 

A.     APPROACH 

The  transfer  function  of  a  two-dimensional  FIR  filter  with  impulse  response 
support  in  the  first  quadrant  can  be  defined  as 

H(zliZ2)=^2     £  h{ni,n2)z;niz-n\  (2.1) 

n  i  =  0    n  3  =  0 

The  frequency  response  of  a  such  a  filter  is  obtained  by  evaluating  the  transfer 
function  on  the  unit  bicircle  through  the  substitutions  zx  —  eJWl  and  z2  =  e3"2 , 
yielding 

ATi-l  Na-1 

F(«Xlwi)=   £)    £  M%>"2)e-ywinie-iwan'.  (2.2) 

The  frequency-sampling  filter  design  problem  involves  determination  of  a  sequence 
h{nl ,  n2 )  which  will  result  in  a  filter  matching  a  desired  frequency  response  Hd  {uj1  ,  w2 ) 
at  specified  sample  points  in  the  (w1?u;2)  plane. 

The  initial  step  in  the  design  process  is  to  choose  the  filter  order,  and  hence 
the  parameters  N^   and  JV2.    Often  it  is  desirable  to  specify  Nt   =  N2   —  N  to 
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ensure  a  result  with  equal  dependence  on  the  frequencies  ujx  and  u>2  •  The  desired 
frequency  response  iyd(wi,w2)  is  sampled  at  the  vertices  of  a  uniformly  spaced 
jVj  x  AT2  Cartesian  grid  on  the  region  {0  <  u^  <  27r;  0  <  u2  <  2n}  of  the  (u^  ,u>2) 
plane.  The  resulting  frequency  samples  can  be  expressed  in  2-D  sequence  form  as 


Arx  =  0, 1, ...  ,7V!  —  1     and     k2  =  0, 1, . . .  yN2  -  1. 
The  resulting  sampling  geometry  is  illustrated  in  Figure  2.1. 


2  rr  f  Af   t   -    1 
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2»r  I  N   i  -    1  ) 


(2.3) 


Figure  2.1.  Arrangement  of  Samples  for  Uniform  Sampling  Method. 


Because  an  IDFT  approach  is  to  be  used  in  computing  h{nx ,  n2),  and  since  the 
IDFT  is  defined  for  first-quadrant  sequences  only,  the  desired  frequency  response 
should  possess  a  linear  phase  characteristic  such  that 


ffd(wi,«2)  =  \Hd(u>l,u>2)\e->ulM*e 


J  W  i  A/  j        -  J  U  j  Af  ; 


(2.4) 
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where 


N,  /2  Nt  even  ,      , 

Mi:     M^-l)/2     A^odd  <2-5> 


and 

f       AT2/2 

\  (iV2  -  l)/2      JV2  odd 

Additionally,  Hd(u>i,ui2)  must  possess  the  symmetry  conditions 


,,  jv2/2  N2  even  ,      . 

M*  =  \  /at        iwo      *r    ^     •  (2-6) 


W,fca)  =|H(*i,iVa  -*3)| 
=  |fT(JVi  —  Arr,/c2)| 


(2.7) 


in  order  to  ensure  a  purely  real  impulse  response. 

The  next  step  is  to  compute  h{nx,n2)  by  taking  the  IDFT  of  the  sampled 
desired  frequency  response: 

h{ni,n2)  =JD¥T[H(kl9k2)] 

w,-i  n7-i  ^2  8) 

1        2    kx-0    kj  =  0 

where 


3  «■ 


WNl  =  e   3  "i      and     W^,  =  e   J  w*  . 

The  resulting  sequence  /i^ ,  n2)  is  purely  real  and  represents  the  impulse  response 
of  the  designed  filter.  Such  a  filter  possesses  a  linear  phase  characteristic.  A  zero- 
phase  characteristic  can  be  obtained  by  shifting  the  impulse  response  so  that  it  is 
centered  about  the  origin  such  that 

h'(mx,m2)  =  h{mx  +  Mx,m2  +  M2),  (2.9) 

where  h'(mx,m2)  represents  the  zero-phase  impulse  response. 
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To  see  the  relationship  between  the  resulting  filter  frequency  response  and  the 
specified  sample  sequence  H(ki,k2),  Eq.  (2.8)  is  substituted  into  Eq.  (2.2): 

m  =  0    n3  =  0   1-12    fci  =  0    kl-0 


JVtAT2 


fc,  =0    fc3  =  0 


n  i=0 


n,  =  0 


JV!-  1  JV3-  1 
&!  =0    fc3  =  0 


(1  -  e"J'WlJVl)(l  -  e_JW'"') 


H{kx,k2) 


NXN2 


(1  -  W^-ifcle-^')(l  -W"  *'«-*«>) 


(2.10) 


The  above  expression  can  be  rewritten  in  the  form 

«,  -  1  W3 -  1 

H(k,.k„)$(u,  - 

Kj=U      Sj  =  U 

where 


#(wl5w2)  =     V     V   #(*!,*;,)$(«!  -  T^*i,wa  -~*a).  (2.H) 

iVi  iV2 

fc1=0    fc3=0 


$(<*>i,w2) 


(1  -  e~JWlNl)(l  -  e_yw5Wj) 
NlN2{\  -  e--""i)(l  -  e-'w») 

•  e-7"1 '     » — '  e~',W3^     a — '  e 


-y«i(i£Vi)e-y«»(iIVli)ei«*i(i-j#r)cy**»(i-j^) 


which  simplifies  to 


$(wl5w2)  = 


sin(w1iV1/2)  sin(u;2iV2/2) 


ATj  JV2  sin(w1  /2)  sin(u;2  /2) 


(2.12) 


g-jwi(    , — )g-y«M         > g- 


for  the  case  where  Nx  and  N2  are  both  odd  integers.  This  analysis  is  a  2-D 
extension  of  operations  performed  for  the  one-dimensional  case  [Ref.  1:  pp.  97- 
98].  Since 


&(—k     —k  )       I  1 
[NX    l'N2    2>  "  1  0     ifc1=l,2J...,iVi--l     and     k2  ==  1,2,. . .  tNa  -  1 ' 


(*i,Aa)  =  (0,0) 

—  1     and     /c2  = 

the  frequency  response  is  matched  exactly  at  the  original  sample  points.  That  is, 


#(<*>!, w2)L  _i^x  M  _i-*i  =  H{kx,k2), 


(2.13) 
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kx  =  0,1,...,^  -  1     and     k2  =  0, 1, . . . ,  N2  -  1. 

Thus,  the  design  process  entails  a  bivariate  interpolation  in  the  variables  u)x  and 
u2,  and,  due  to  use  of  the  IDFT  algorithm,  the  interpolating  functions  are  of  the 
form  of  Eq.   (2.12). 

In  using  this  filter  design  method,  the  choice  of  the  desired  frequency  response 
Hd[ujl  ,u>2)  has  significant  effects  on  the  quality  of  the  resulting  filter.  For  instance, 
if  an  "ideal"  lowpass  filter  characteristic  with  the  circular  symmetry  shown  in 
Figure  2.2  is  used  for  Hd  (wj ,  u>2 ) ,  the  passband  and  stopband  of  the  resulting  filter 
frequency  response  exhibit  substantial  ripple. 


Figure  2.2.  Ideal  Lowpass  Filter  Frequency  Response. 
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Figure  2.3  illustrates  contour  and  perspective  plots  of  the  frequency  response 
arising  from  such  a  design.  In  this  case,  a  15  x  15  point  filter  was  designed  by 
uniformly  sampling  the  desired  frequency  response  of  Figure  2.2. 

The  frequency  response  of  the  resulting  filter  can  be  improved  considerably  if 
a  transition  region  is  specified  for  Hd(u)i  ,  u>2).  Increasing  the  width  of  such  a  tran- 
sition region  reduces  passband  and  stopband  ripple  but  does  so  at  the  expense  of  a 
less  sharp  roll-off  characteristic  in  the  resulting  filter  response.  The  values  assigned 
to  H{ki,k2)  for  points  (ki,k2)  located  in  this  transition  region  significantly  influ- 
ence the  behavior  of  the  resulting  filter.  The  selection  of  these  values  in  order  to 
satisfy  some  specific  optimization  criterion  for  the  frequency  response,  such  as  best 
equi-ripple  approximation  to  an  ideal  response,  is  not  trivial.  Since  H(oj1,uj2)  is  a 
linear  function  of  the  values  H[kl^k2),  as  seen  in  Eq.  (2.10),  linear  optimization 
techniques  can  be  employed  in  selecting  the  transition  region  values  of  H(ki,h2). 
Such  techniques  for  the  one-dimensional  case  are  rather  involved  [Ref.  4];  for  the 
2-D  case,  the  problem  is  much  more  computationally  intensive  [Ref.  2],  and  hence, 
has  not  been  addressed.  Specifying  Hd(u>i,<jj2)  to  have  a  smooth,  linear  transition 
between  the  passband  and  stopband  regions  is  the  simplest  thing  to  do  and  gen- 
erally works  well.  The  magnitude  characteristic  of  a  desired  frequency  response 
with  a  linear,  circularly  symmetric  transition  region  is  shown  in  Figure  2.4.  A  fil- 
ter can  then  be  designed  by  sampling  this  response  on  the  uniform  Cartesian  grid 
specified  in  Eq.  (2.3).  Contour  and  perspective  plots  of  a  15  x  15  filter  designed 
in  this  manner  are  shown  in  Figure  2.5.  In  comparing  this  result  to  that  of  Figure 
2.3,  it  is  evident  that  inclusion  of  a  transition  region  greatly  reduces  passband  and 
stopband  ripple. 
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Figure  2.3.  Frequency  Response  Arising  From  Sampling  Ideal  Lowpass 
Characteristic  With  No  Transition  Band. 
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Figure  2.4.  Desired  Lowpass  Frequency  Response  Characteristic 
With  Linear  Transition  Band. 


B.      EXISTENCE  AND  UNIQUENESS 

A  significant  property  of  the  uniform  sampling  design  is  that  a  unique  filter 
satisfying  each  of  the  interpolating  conditions  specified  in  Eq.  (2.3)  is  guaranteed 
to  exist.  In  general,  existence  and  uniqueness  of  such  an  interpolation  result  in  two 
dimensions  is  not  guaranteed;  this  can  be  a  significant  theoretical  and  practical 
problem.  The  point  is  discussed  in  further  detail  in  Chapter  III.  A  proof  of  existence 
of  a  unique  h(n1,n2)  using  the  frequency-sampling  method  with  uniformly  spaced 
samples  follows. 
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Figure  2.5.  Frequency  Response  Arising  From  Sampling  Desired  Lowpass 
Characteristic  With  Linear  Slope  Transition  Band. 
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Equation  (2.8)  is  used  to  determine  h(n1,n2).  This  can  be  rewritten  as 


1     fc1=0L       2    fc3  =  0 


fc-jn- 


(2.14) 


Now  define 


jv,  - 1 


G{kl,na)  =  jr   J2  H(kl:k2)W-^ 


(2.15) 


fc,  =  0 


Substituting  Eq.    (2.15)  into  Eq.    (2.14),  the  following  expression  for  the  impulse 
response  results: 

M»i.«a)  =  T7-   £  G^x.najW-*1"1.  (2.16) 


#i 


fc1  =  0 


Thus  the  2-D  IDFT  of  H(k1 ,  k2)  can  be  determined  in  a  two-step  process  involving 
one-dimensional  transforms  known  as  row-column  decomposition  [Ref.  5:  pp.  75- 
76]. 

The  first  step  in  the  row-column  decomposition  process  is  to  perform  the 
solution  of  a  total  of  Nx  systems  of  N2  equations  of  the  form  of  Eq.  (2.15)  using  a 
1-D  IDFT  algorithm.  In  matrix  form,  this  can  be  written  as 


G(*i,0) 
G(*i,l) 


N7 

1 

w^l 

iv7 

• 

• 

w 


N2 
-(Ni-1) 


N- 


lG{kx,N2-l) 

which  can  be  expressed  in  vector-matrix  notation  as 


H(kly0) 
H(*i,l) 


W-1      J    IHik^Ni-l)] 

(2.17) 


G  =  -W2H. 


(2.18) 


It  can  be  shown  [Ref.  6:  pp.  44-45]  that  a  scaled  version  of  W2  ,  namely  W2  /y/N2i 
is  unitary  and  hence  non-singular.  Therefore  a  unique  solution  for  G  exists  for  each 
k,      t  =  0,l,...,#x  -I- 
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The  second  step  in  the  computation  involves  the  solution  of  7V2  systems  of  Nt 
equations  of  the  form  of  Eq.  (2.16).  In  matrix  form,  this  becomes 


h(0,n2) 
h{l,n2) 

lh(Nx  -l,n2)J 


AT, 

1 

w° 

N~i 

! 

\ 

LKKJV-1 

w 


Nx 
-[Nx-X] 


AT, 


G{0,n2) 
G(l,n2) 


W-1       J    IGiN,  -l,n2)J 

(2.19) 


which  in  vector-matrix  form  is 


h=-WlG. 


(2.20) 


Since  G  was  uniquely  determined  for  each  value  of  kXi  and  since  Wj  is  also 
non-singular,  a  unique  solution  for  h  exists  for  each  ny,  j  =  0, 1, ... ,  JVjj  —  1. 
Therefore,  a  unique  solution  for  h(nl,n2)  is  obtained  from  the  sequence  H(kl,k2), 
which  represents  the  samples  of  a  desired  frequency  response  at  the  vertices  of  a 
uniform  Cartesian  grid  in  the  (wl5w2)  plane. 

C.      COMPUTATIONS 

While  there  is  no  flexibility  in  the  choice  of  sample  point  locations  using  this 
IDFT  approach,  the  method  does  offer  some  significant  computational  advantages. 
One  criterion  often  used  to  measure  computational  efficiency  in  solving  a  problem 
is  to  determine  the  number  of  complex  operations  involved.  While  there  exist  many 
methods  of  computing  the  2-D  IDFT,  the  use  of  row-column  decomposition  and 
the  IDFT  directly  (vice  the  FFT  algorithm)  is  a  reasonably  efficient  approach  for 
the  design  of  filters  of  practical  order  and  is  used  in  the  following  analysis. 

The  initial  step  in  the  row-column  decomposition  approach  is  to  solve  a  total 
of  Nx  systems  of  the  form  of  Eq.  (2.17).  Solution  of  each  system  via  the  IDFT 
involves  a  total  of  JV2   multiplications  and  additions.   Therefore,  a  total  of  NiN2 
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floating  point  operations  are  required  in  this  step,  where  a  floating  point  operation 
is  considered  one  multiplication  and  one  addition. 

The  second  part  of  the  row-column  decomposition  approach  involves  the  so- 
lution of  N2  systems  of  the  form  of  Eq.  (2.19).  Solution  of  each  system  via  the 
IDFT  requires  JVf  operations.  Therefore,  Nf  N2  floating  operations  are  involved 
in  the  second  step. 

The  total  number  of  floating  point  operations  required  in  the  solution  for  the 
impulse  response  coefficients  of  an  Ni  x  N2  frequency-sampling  FIR  filter  design 
using  uniform  sampling  on  a  Cartesian  grid  and  a  row-column  decomposition,  IDFT 
approach,  is 

Cuniform  =  N1N^+N^N2.  (2.21) 

For  the  case  in  which  Nx  —  N2  =  N,  this  reduces  to 

Cuniform  =27V\  (2.22) 

In  summary,  the  uniform  sampling,  IDFT  approach  to  2-D  frequency  sampling 
filter  design  has  the  benefit  of  computational  efficiency,  yet  lacks  flexibility  in  the 
choice  of  frequency  samples.  The  next  two  chapters  will  present  methods  which 
achieve  greater  flexibility  by  allowing  nonuniform  spacing  of  samples. 
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III.     DESIGNS  WITH  ARBITRARILY  PLACED  SAMPLES 

In  Chapter  II  a  frequency-sampling  2-D  FIR  filter  design  method  was  pre- 
sented in  which  samples  were  placed  on  a  uniform  Cartesian  grid  in  the  (oj1,oj2) 
plane.  Due  to  the  constraints  on  sample  placement,  however,  the  method  was  not 
particularly  flexible.  If  the  desired  frequency  response  Hd  (u^  ,  u>2)  was  fixed,  no  free 
parameters  were  available  in  optimizing  the  frequency  response  of  the  resulting  fil- 
ter. In  this  chapter,  an  alternate  frequency-sampling  design  technique  for  2-D  FIR 
filters  is  presented  in  which  samples  are  no  longer  constrained  to  a  uniform  Carte- 
sian grid  but  can  be  placed  anywhere  in  the  (oj1,oj2)  plane.  This  approach  will 
be  referred  to  as  the  "arbitrary-sampling"  method.  Allowing  the  samples  to  occur 
anywhere  in  the  (wx,c«;2)  plane  introduces  additional  free  parameters;  hence,  linear 
optimization  techniques  can  be  employed  without  varying  Hd  (wj ,  u2 ) ,  the  desired 
frequency  response  to  be  sampled.  Since  the  locations  of  the  frequency  samples  are 
much  less  constrained,  the  method  provides  greater  potential  for  achieving  a  desir- 
able frequency  response  characteristic  in  the  resulting  filter.  Since  the  placement 
of  samples  on  a  uniform  Cartesian  grid  is  a  special  case  of  the  arbitrary-sampling 
method,  the  arbitrary-sampling  method  can  achieve  results  at  least  as  good  as 
the  uniform-sampling  method  for  a  fixed  impulse  response  region  of  support  and 
a  fixed  number  of  samples.  However,  the  arbitrary-sampling  method  possesses  a 
few  significant  drawbacks,  such  as  computational  complexity  and  the  possibility  of 
degenerate  cases  in  which  the  interpolating  conditions  cannot  be  met.  These  must 
be  considered  in  the  selection  of  a  method  for  filter  design. 
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A.     APPROACH 

The  arbitrary-sampling  method  can  be  used  to  produce  2-D  FIR  filters  with 
a  variety  of  impulse  response  regions  of  support.  This  region  of  support  is  one  of 
the  parameters  to  be  specified  at  the  outset  of  the  design  process.  For  illustration 
purposes,  a  filter  with  an  Nx  x  7V2  rectangularly-shaped  impulse  response  region 
of  support  centered  on  the  origin  of  the  {nl,n2)  plane  is  specified,  where  Nt  and 
N2  are  odd  integers.  The  transfer  function  of  such  a  filter  can  be  expressed  as 

H{z1,z2)=      J2  L      H^,n2)z;n^z-n\  (3.1) 

where  h(n1,n2)  is  the  impulse  response  to  be  determined,  Mx  =  (Ni  —  l)/2  and 
M2  =  (N2  —  l)/2.  The  parameters  Mi  and  M2  will  be  used  often  in  the  remainder 
of  this  analysis,  so  the  above  relationships  are  worth  noting.  Evaluation  of  Eq. 
(3.1)  on  the  unit  bicircle  yields  the  following  expression  for  the  frequency  response: 

Mx  M, 

J5T(wi,wa)=      £         Y,     M^i^2)e-iWxnie-JWan3.  (3.2) 

n  i  =  —  M  i  ti)  =  —  M  -j 

To  ensure  a  zero-phase  characteristic  in  the  resulting  frequency  response, 
h{n1,n2)  is  specified  to  possess  quadrant  symmetry  such  that 

/i(n!,n2)  =  h(—nltn2)  =  h(n1,—n2).  (3.3) 

Application  of  Eq.  (3.3)  to  Eq.  (3.2)  produces  the  following  expression: 

Mx 

-H"(u>i,w2)  =M°>°)  +  Yl  2Mnn°)coswini  + 


n  x  =  1 

hi  2  hi  i  hi  ; 


V^  2h(0,n2)  cos u>2n2  +   YJ    YJ  Ah{nx  ,n2)  cos ^rij  cosw2n2. 

(3.4) 


fi2  =  l  Ri  =  l  ni  =  l 
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This  can  be  rewritten  as 

M1  M; 


H((jj1,u>2)  =    /_.    /v  a(^i,w2)cosw1n1cosu;2n2,  (3.5) 


n  i  =  0  ni  =  0 


where 

Ih(n1,n2)  n1  =  0,n2  =  0 

2/i(nl5n2)  n!^0,ri2=0  r      * 

2Mnlin9)  n1=0,n2^0    '  (3"6) 

4/i(ni ,  n2 )  ^^0,^2^0 

The  basis  functions  {t/;,  (c^!  ,w2)}  are  defined  such  that 

■0,  (wi  ,u;2)  =  cos^n!  coscj2ra2,  (3-7) 

where 


(3.8) 


i  =  nl  (M2  +  1)  +  n2  +  1;         r^  =  0, 1, . . .  ,MX ; 

n2  =  0, 1,...,M2. 
By  arranging  the  2-D  sequence  a{nx,n2)  into  the  one-dimensional  sequence  a(i), 

the  frequency  response  can  then  be  expressed  as 

R 

ff(wi,w2)  =  ^a(#(w1,w2),  (3.9) 

i=  1 

where 

R  =  {MX  +1)(M2  +1).  (3.10) 

From  Eq.  (3.9),  the  elements  of  the  sequence  a{i)  represent  a  total  of  R  free 
parameters  to  be  determined.  Therefore,  the  frequency  response  is  to  be  specified 
at  a  total  of  R  distinct  sample  points  in  the  (oj1,u>2)  plane. 

The  initial  step  in  the  design  process  is  to  specify  H[ui,uj2)  at  R  distinct 
sample  points  in  the  region  K  —  {{oJ\,oj2)  :  0  <  ux  <  7r;0  <  u2  <  tt}.  An  example 
of  such  a  sampling  geometry  is  illustrated  in  Figure  3.1. 

No  constraints  other  than  the  distinctness  requirement  are  placed  on  the  loca- 
tion of  samples  within  this  region.  The  form  of  Eq.  (3.5)  implies  that  the  frequency 
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Figure  3.1.  Possible  Arrangement  of  Samples  for  Arbitrary  Sampling  Design. 

response  possesses  quadrantal  symmetry  and  periodicity  of  2n  in  both  frequencies. 
Hence,  /f(wi,w2)  is  defined  for  all  (u^Wj). 

The  next  step  in  the  design  process  is  to  solve  for  the  sequence  a(ni  ,n2).  Appli- 
cation of  Eq.  (3.9)  at  each  specified  sample  location  {(wj  ,u;2),  ;  i  =  1,2, . . .  ,/2} 
results  in  the  following  system  of  linear  equations: 


#(wi,w2)i  " 
#(<*>i  ,w2)2 

H(u{,u2)r  - 


^1(^1,^2)1  02(Wi,W2)i 

^1(^1,^2)2    ?Mwi,w2)3 


V>«(wi,w2)i 

V»*  (^1,W2)2 


-V,i(wi,w2)fl      02(^1,^2)/?      •••     ^(^1,^2)* 


a(l) 

a(2) 


a(R) 
(3.11) 
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In  vector-matrix  notation,  Eq.  (3.11)  becomes 

H  =  Ba,  (3.12) 

in  which  B  is  of  dimension  R  x  R.  The  sequence  0,(1%!  ,n2)  results  from  solving  the 
relation 

a  =  B-xH.  (3.13) 

The  final  step  in  the  process  is  to  compute  the  filter  impulse  response  h{nl ,  n2) 

using  the  relationship  of  Eq.  (3.6).  The  resulting  impulse  response  has  an  Nx  x  N2 

region  of  support.  Equation  (3.4)  suggests  a  transfer  function  of  the  form 

M  x  m  , 

H{zl,z2)=h{0,0)+  Y,  Mwx.OJW  +z;ni)+  J2  M0,n2)(*2S  +^"3) 

n  i  =  1  nj  =  l 

Mi        Mi 

(3.14) 
which  can  be  used  in  implementing  a  2-D  filter  structure  with  a  minimal  number 
of  gains. 

This  same  approach  can  also  be  employed  if  other  than  rectangular  impulse 
response  regions  of  support  are  specified.  Again,  a  system  of  linear  equations  of 
the  form  of  Eq.  (3.5)  is  solved,  but  the  limits  on  the  summations  are  altered  to 
reflect  the  desired  region  of  support  for  h(rii  ,n2). 

To  summarize  the  arbitrary-sampling  design  process,  the  following  steps  are 
performed: 

(1)  Fix  the  impulse  response  region  of  support,  and  hence,  the  filter  order. 

(2)  Specify  the  desired  frequency  response  at  selected  sample  points  on  the 
region  K  —  {(u^ ,  w2)  :  0  <  ujx  <  7r;  0  <  w2  <  7r}.  The  number  of  specified 
samples  must  equal  the  number  of  independent  impulse  response  values 
(i.e.,    those  which  are  not   fixed  by   the  symmetry  conditions).      The 
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locations    of   samples    within    the    region    is    unconstrained    except    for 
the  condition  that  each  sample  location  be  distinct. 

(3)    Solve    a   system   of   linear   equations   to   obtain   the    resulting   impulse 
response  h{nl ,  n2 ) . 

B.     EXISTENCE  AND  UNIQUENESS 

There  is  a  fundamental  deficiency  in  the  arbitrary-sampling  approach  to  the  2- 
D  FIR  frequency-sampling  filter  design  problem.  For  specified  frequency  response 
samples  placed  arbitrarily  in  region  K,  there  is  no  guarantee  that  a  unique  filter 
possessing  a  number  of  independent  impulse  response  samples  equal  to  the  number 
of  interpolating  conditions,  and  which  satisfies  each  of  those  interpolating  condi- 
tions, exists.  This  is  a  problem  inherent  to  interpolation  in  two  dimensions  which 
does  not  arise  in  the  one-dimensional  case. 

For  a  unique  solution  to  the  system  of  equations  indicated  by  Eq.  (3.12)  to 
exist,  the  basis  functions  {ipt :(w1  ,  cu2)}  must  satisfy  the  Haar  condition.  A  set  of 
vectors  in  n-space  satisfies  the  Haar  condition  if  every  set  of  n  of  them  is  linearly 
independent  [Ref.  7:  p.  45].  Therefore,  to  satisfy  the  Haar  condition,  any  set  of 
R  characteristic  vectors  defined  on  R  distinct  points  in  region  K  of  the  (c^,^) 
plane  must  be  linearly  independent.  For  the  arbitrary-sampling  method  developed 
in  this  chapter,  a  characteristic  vector  {$((wi,w2)i)}  associated  with  a  particular 
sample  location  (w^Wj);  can  be  defined  as 

*((wi,Wa)*)  =  [lM(wi»wa)«)      W(wi,u>2)i)      •••     ^R  ((wi  »  wa)i)  f  •        (3-15) 

It  is  known  that  for  a  region  K'  of  the  two-dimensional  frequency  plane,  there  are  no 
nontrivial  sets  of  two-dimensional  basis  functions  which  satisfy  the  Haar  condition 
[Ref.  8:  p.  494].  Therefore,  there  is  no  guarantee  that  the  sets  of  characteristic 
vectors  formed  by  evaluating  Eq.    (3.15)  at  arbitrarily  selected  sample  points  in 
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region  K  of  the  (uhw2)  plane  are  linearly  independent.   Since  the  square  matrix 
B  in  Eq.   (3.12)  can  be  expressed  in  terms  of  the  characteristic  vectors  as 


B  =  [tBr((w1,w3)1)     9{{ultu2)a)     ...     *((w1)w2)fl)]r, 


(3.16) 


and  since  the  characteristic  vectors  may  be  linearly  dependent,  B  may  be  singular 
and  a  solution  to  Eq.  (3.13)  may  not  exist.  Such  a  case  is  termed  a  degeneracy, 
and  illustrates  a  basic  problem  associated  with  the  arbitrary-sampling  method.  As 
can  be  inferred  from  the  form  of  the  basis  functions,  the  occurrence  of  a  degeneracy 
is  dependent  solely  upon  the  location  of  samples;  the  specified  frequency  response 
values  at  the  sample  points  do  not  affect  this. 

A  simple  example  illustrates  how  a  degeneracy  can  arise.  The  frequency  re- 
sponse of  a  filter  is  specified  at  four  sample  points  in  the  (u!x  ,u2)  plane: 

1.  H(0, 0.6tt)  =  1  2.  Zr(0.47r,0)  =  1 

3.  #(0.47r,7r)  =  0  4.  i7(7r,0.67r)  =0 

These  sample  locations  are  shown  in  Figure  3.2. 

The  parameters  Mx  and  M2  are  both  set  equal  to  1  so  that  the  number  of 
samples  in  the  sequence  a{n1,n2)  equals  the  number  of  frequency  samples.  Sub- 
stitution into  Eq.  (3.11)  yields 


1 

1 

0 

LOJ 


1       COS0.67T  1  COS0.67T 

1  1  COS0.47T  COS0.47T 

1  —1  COS0.47T       —  COS0.47T 


a(0,0) 
a(0, 1) 
a(l,0) 


(3.17) 


.1     cos0.67r  —1  —  cos0.67rj    La(l,l). 

Some  calculations  reveal  that  the  square  matrix  above  has  a  determinant  of  zero; 
therefore,  it  is  not  invertible  and  no  unique  solution  for  a{rii,n2)  exists. 

To  determine  whether  degenerate  cases  present  a  serious  shortcoming  of  this 
method,  it  is  of  interest  to  determine  how  often  degeneracies  are  likely  to  arise.  To 
investigate  this,  an  APL  computer  program  was  developed  to  count  the  number  of 
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Figure  3.2.  Arrangement  of  Samples  for  Design  Example. 

degenerate  cases  arising  from  the  arbitrary  selection  of  sample  points  {(c*>i  ,w2)i} 
on  region  K  of  the  (wi,w2)  plane.  A  total  of  R  samples  were  arbitrarily  placed 
at  distinctly  different  vertices  of  a  grid  of  specified  density  on  region  K  for  each 
case.  The  placement  of  samples  was  performed  in  a  pseudo-random  fashion  using 
the  APL  "roll"  function.  The  density  of  the  grid,  in  effect,  is  determined  by  the 
number  of  decimal  places  carried  in  the  frequency  values.  In  the  analysis,  two 
parameters  were  varied:  the  dimensions  of  the  specified  filter  impulse  response, 
and  the  grid  density.  For  each  situation,  one  hundred  separate  sets  of  arbitrary 
samples  were  investigated.  The  program  then  determined  the  number  of  cases  in 
which  the  sampling  geometry  resulted  in  a  degeneracy.  Due  to  the  large  number  of 
computations  involved,  analysis  was  limited  to  filters  with  square  impulse  response 
regions  of  support  only.  The  results  are  listed  in  Table  3-1. 
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These  results  support  two  concepts.  First,  for  a  fixed  grid  density,  as  the 
order  of  the  specified  filter  increases  (and  hence  the  number  of  specified  frequency 
samples),  the  likelihood  of  the  occurrence  of  a  degeneracy  increases.  Secondly,  for  a 
fixed  filter  order,  as  the  density  of  the  grid  on  which  samples  are  placed  increases, 
the  likelihood  of  a  degeneracy  occurring  decreases.  For  filters  of  order  17  x  17 
or  less,  the  probability  of  a  degenerate  case  is  minimal  provided  frequencies  are 
specified  to  at  least  the  nearest  0.05  radians. 

TABLE  3-1 
NO.  OF  DEGENERACIES  OCCURRING  IN  100  TEST  CASES 


GRID  DENSITY  (RADIANS) 

FILTER  ORDER 

R 

0.2                      0.1                      0.05 

5x5 

9 

5 

1 

0 

9x9 

25 

10 

0 

0 

13  x  13 

49 

17 

1 

0 

17  x  17 

81 

56 

2 

0 

It  is  significant  to  note  that  for  the  above  tests,  frequency  samples  were  cho- 
sen in  an  arbitrary  manner.  However,  when  designing  a  filter  to  approximate 
a  certain  desirable  frequency  response  characteristic,  the  selection  of  points  is 
likely  to  be  done  in  an  unconstrained,  but  not  arbitrary,  manner.  Certain  pat- 
terns in  sample  selection,  such  as  selecting  octant  symmetric  sample  points  where 
(w1,w2)t  =  {w2,0l>i)j,  have  a  greater  tendency  to  introduce  linear  dependencies 
within  the  set  of  characteristic  vectors  than  when  sample  points  are  selected  at 
random.  This  concept  is  difficult  to  show  in  a  controlled  experiment;  however,  in 
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the  course  of  this  research,  degeneracies  arose  when  attempting  filter  designs  much 
more  often  than  the  arbitrarily  placed  sample  test  results  would  indicate. 

Such  a  degeneracy  did  arise  in  the  following  design  example.    In  this  case,  a 
lowpass  filter  is  desired  with  a  circularly  symmetric  frequency  response  such  that 


it  possesses  unity  gain  in  the  passband  region  {0  <  \/u>*  +  u\   <  0.27r}  and  zero 


gain  in  the  stopband  region  {0.47T  <  \/uj]  +  u\  <  7r}.  The  region  of  support  for 
the  impulse  response  is  specified  to  be  17  x  17,  implying  Mx  —  M2  =  8.  Therefore, 
a  total  of  R  —  [Mx  -f  1)(M2  +1)  =81  frequency  samples  are  to  be  specified.  Figure 
3.3  shows  the  location  of  81  points  at  which  the  desired  frequency  response  is  to 
be  sampled  initially. 


-o  u 


OJ 


Figure  3.3.  Arrangement  of  Samples  for  Lowpass  Filter 
Design,  Case  1. 
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Note  that  the  arbitrary-sampling  method  does  not  require  any  samples  to  be 
placed  in  the  transition  (don't  care)  region.  This  allows  the  concentration  of  more 
samples  in  the  passband  and  stopband  regions.  This  particular  choice  of  sample 
locations,  however,  was  found  to  result  in  a  degeneracy  when  a  computer  solution 
was  attempted. 

Following  this  degenerate  result,  the  sampling  geometry  was  altered  by  slight 
displacement  of  the  locations  of  six  of  the  frequency  samples.  The  new  sampling 
geometry  is  illustrated  in  Figure  3.4,  with  the  circled  points  indicating  the  sample 
locations  which  differ  from  the  degenerate  case.  The  resulting  filter  frequency 
response  is  shown  in  Figure  3.5. 


L 


10 


D  : 


Figure  3.4.  Arrangement  of  Samples  for  Lowpass  Filter 
Design,  Case  2. 
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Figure  3.5.  Frequency  Response  of  Lowpass  Filter 
Design,  Case  2. 

Although  the  result  may  seem  to  be  perfectly  acceptable,  the  system  of  equa- 
tions involved  in  the  solution  can  be  considered  to  be  ill-conditioned.  Slight  per- 
turbation of  the  position  of  a  single  sample  point,  as  shown  by  the  circled  point  in 
Figure  3.6,  results  in  the  unacceptable,  excessively  rippled  response  indicated  by 
Figure  3.7. 

The  above  example  demonstrates  a  negative  aspect  of  the  arbitrary-sampling 
approach.  Although  degeneracies,  as  defined,  may  seldom  occur,  ill-conditioned 
cases  in  which  the  determinant  of  matrix  B  in  Eq.  (3.12)  is  very  small  ,  but  non- 
zero, may  arise  much  more  frequently.  While  specifying  the  frequencies  u^  and  u2 
at  each  sample  point  with  a  high  degree  of  precision  (analogous  to  placing  samples 
on  a  very  dense  grid)  may  assure  a  very  minimal  probability  of  a  degeneracy,  it 
does  not  guarantee  that  the  resulting  system  will  not  be  ill-conditioned. 
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Figure  3.6.  Arrangement  of  Samples  for  Lowpass  Filter 
Design,  Case  3. 

C.     COMPUTATIONS 

A  further  negative  aspect  of  this  arbitrary-sampling  method  is  the  number 
of  computations  involved  in  determining  filter  parameters.  While  the  uniform- 
sampling  method  of  Chapter  II  was  shown  to  be  reasonably  efficient  in  computation, 
the  arbitrary-sampling  method  involves  the  solution  of  a  rather  large  system  of 
linear  equations.  This  solution  requires  the  inversion  of  an  R  x  R  matrix.  This 
inversion  dominates  the  computation  process.  It  is  well  known  that  inverting  an 
N xN  matrix  requires  on  the  order  of  iV3  floating  point  operations,  where  a  floating 
point  operation  consists  of  one  multiplication  and  one  addition.  Since,  for  a  filter 
with  an  Nt  x  N2  impulse  response  region  of  support, 


R  =  {M1  +  1)(M,  +  1)  =  (2L-1  +  l)(^5-i  +  1) 

=  i(JV,  +  1)(JV,  +  1) 

4 


(3.18) 
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Figure  3.7.  Frequency  Response  of  Lowpass  Filter 
Design,  Case  3. 

the  order  of  the  number  of  floating  point  operations  involved  in  the  arbitrary- 
sampling  design  process  can  be  represented  as 


1 
64 


Carb  *  f^VjJVj. 


(3.19) 


For  the  case  in  which  Nx  =  N2  =  N,  this  reduces  to  the  following: 


Carb  ~  —  N6, 
64 


(3.20) 


Thus,  the  arbitrary  sampling  method  represents  an  approach  which,  while 
quite  flexible,  is  computationally  intensive  and  fails  to  provide  guarantee  of  a 
unique  solution.  In  the  next  chapter,  an  alternative  nonuniform  sampling  approach 
is  presented  which,  by  placing  some  mild  constraints  on  frequency  sample  location, 
always  provides  a  unique  design  solution  and  is  much  more  computationally  effi- 
cient. 
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IV.     A  NOVEL  2-D  NONUNIFORM  FREQUENCY-SAMPLING  APPROACH 

Chapters  II  and  III  presented  two  fundamentally  different  approaches  to  2-D 
FIR  frequency-sampling  filter  design.  The  uniform-sampling  method  possessed  no 
flexibility  in  the  location  of  frequency  samples,  yet  was  computationally  efficient  in 
determining  filter  parameters.  The  arbitrary  sampling  method,  however,  had  few 
constraints  on  sample  location  but  was  computationally  intensive.  The  approach 
described  in  this  chapter  represents  a  compromise.  It  is  more  constrained  than  the 
arbitrary  sampling  method,  but  has  much  more  freedom  in  sample  placement  than 
for  the  uniform-sampling  method.  While  not  as  efficient  as  the  IDFT  approach, 
this  method  requires  significantly  fewer  computations  when  compared  with  the 
arbitrary  sampling  method.  Additionally,  existence  and  uniqueness  of  a  solution 
for  the  impulse  response  is  guaranteed,  provided  simple  constraints  on  sample 
locations  are  met. 

A.     APPROACH 

For  convenience,  the  design  method  introduced  in  this  chapter  will  be  referred 
to  as  the  "nonuniform  sampling"  method,  although  the  arbitrary  sampling  method 
uses  nonuniformly  spaced  samples  as  well.  This  method  borrows  ideas  from  the 
other  two  approaches.  As  in  the  arbitrary-sampling  case,  quadrant  symmetry  is 
specified.  As  was  shown  in  Chapter  III,  the  frequency  response  of  a  2-D  zero-phase 
FIR  filter  with  an  JVi  x  N2  impulse  response  region  of  support  can  be  represented 
as 

H{fjJx,ui2)  =   2.    /_,  a{ni  >n2)  cosujitii  cosu>2n2,  (4-1) 


ni=0  ni  =  0 
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where 

Ih(nl,n2)  nx  —  Q,n2  =0 

2h(nx,n2)  nx  ^  0,n2  =  0  , 

2Mnx.fi,)  n1=0,n2^0'  ^ 

4/i(n!,n2)  ni^0,n2^0 

with  Mx  =  (Ni  —  l)/2  and  M2  =  (N2  —  l)/2.  In  the  arbitrary  sampling  method, 
the  sequence  a(n1,n2)  was  represented  as  a  column  vector,  the  basis  functions 
^i(w1,w2)  were  specified,  and  a  system  of  equations  of  the  form  of  Eq.  (4.1)  was 
solved  directly.  In  the  nonuniform  sampling  method,  constraints  placed  on  the 
sampling  geometry  result  in  Eq.  (4.1)  becoming  a  separable  system.  The  system 
is  then  solved  by  a  two-step  process  analogous  to  the  row-column  decomposition 
used  in  the  uniform  sampling  method  of  Chapter  II.  The  resulting  systems  of  linear 
equations  contain  one-dimensional  basis  functions  only;  hence,  linear  dependencies 
associated  with  two-dimensional  basis  functions  are  not  of  concern.  This  results  in 
assurance  of  a  unique  solution  to  the  filter  design  problem. 

As  in  the  arbitrary  sampling  method,  samples  are  to  be  placed  in  region  K  = 
{(u>i,u>2)  :  0  <  Wi  <  7r;0  <  uj2  <  7r}  of  the  two-dimensional  frequency  plane. 
As  mentioned  above,  certain  constraints  on  the  location  of  frequency  samples  are 
required.  Suppose  any  (M:  +  1)  distinct  values  of  u1  are  chosen  on  the  interval 
0  <  Wi  <  7r.  These  values  are  denoted  as  {u;lfc  :  k  =  0, 1, . . .  ,MX }.  Then  for  each 
wlk,  suppose  any  (M2  +  1)  distinct  values  of  w2  are  chosen  such  that  0  <  w2  <  7r. 
These  values  are  denoted  as  {uj21  (k)  :  /  =  0, 1, . . . ,  M2  }.  Then  a  unique  solution  for 
a(n1,n2)  can  be  obtained  from 

H{ulk,U2i(k))  =  ^    ^2  a(nlin2)cosulkn1cosu2i(k)n2\ 

0<A;<Ml5     0</<M2. 
Such  a  sampling  geometry  is  shown  in  Figure  4.1  for  the  case  Mx   —  M2   —  2. 
Notice  that  the  selection  of  sample  points  is  quite  flexible.    The  values  w10,  Wn, 
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and  cj12  can  be  any  arbitrary,  distinct  samples  of  u)x  in  the  region.  Additionally, 
w2o(A:)i  w2i(^),  and  uj22{k)  represent  any  arbitrary,  distinct  samples  of  u>2  for  each 
k,  k  —  0,1,2.  This  is  far  less  constrained  than  placing  the  samples  on  a  uniform 
Cartesian  grid. 


u>. 


'*  i 

i 

71 

•  w21(0) 

•  w22(l) 

•  u;22(2) 

•  ^2l(2) 

•    w2i(l) 

• 

•    W20(0) 

•    ^2o(l) 

•  w20(2) 
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0 


1  0 
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1  1 


7T 
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Figure  4.1.  Nonuniform  Sampling  Example. 

Given  the  arrangement  of  samples  described  above,  the  sequence  a(ni  ,  n2)  can 
be  found  as  follows.  First,  Eq.  (4.1)  is  arranged  as 

Mj  Mi 


H{ul,u2)  =   2.     Z.  a(n  1,n2)  cos  uJiUx 

n,=0     n  ,  =  0 

A  new  function  g[u>i,n2)  is  defined  as 

AY, 


cos  ui2  n2 . 


(4.4) 


g{ul,n2)  =   2_.  o.{nl,n2)  cos  ujlnl 


(4.5) 


n,  =0 
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The  frequency  response  can  then  be  expressed  as 

M  j 

ff(wi,w2)  =   2J  ff(wi,nj)cosw2n2. 
For  a  particular  value  of  u>i ,  namely  wlfc ,  the  response  can  be  written  as 

Mi 

H(ulk,u2)  =   ^2  g{uik,n2)cosu2n2. 

ni  =  0 

Now,  for  each  u>2i  (&),  /  =  0, 1, . . . ,  M2 ,  the  above  expression  becomes 
H(ojlk,oj2l{k))  =   2J  g(wifc,n2)cosu;2,(/i:)n2. 


(4.6) 


(4.7) 


(4.8) 


n-,  =  0 


Thus,  the  first  step  in  the  row-column  decomposition  is  to  solve  a  system  of  (M2  +1) 
linear  equations  of  the  form  of  Eq.  (4.8).  In  matrix  form,  this  becomes 


H(ujlk,uj2o{k)) 
H(ulk,uj21(k)) 


1      cosw20(A:)       ...      cosM2w20(A:) 
1       cosu;21(A;)       ...       cosM2o;21  (A;) 


H(u>ik,uj2M2(k))}        Ll     cosui2M:t(k)     .. 

In  vector-matrix  notation,  Eq.  (4.9)  is  expressed  as 

Hfc  =  Agfc , 


g{uik,o) 
0(wlfc,i) 


cosM2w2M](fc)  J    Lg(u>lk,M2). 

(4.9) 


(4.10) 


where  Hfc  and  gfc  are  (M2  +  1)  x  1  column  vectors,  and  A  is  an  (M2  +  l)  x  (M2  + 1) 
square  matrix.  The  vector  gfc  is  then  computed  using  the  relation 


gfc=A-xH 


(4.11) 


This  entire  process  is  performed  for  each  k,  k  =  0, 1, ... ,  M1 .  Therefore,  this  step 
involves  a  total  of  (Mx  +  1)  systems  of  (M2  +  1)  linear  equations  each. 

Once  the  column  vectors  gfc  have  been  computed  for  all  values  of  k,  a  matrix 
G  of  size  (M2  +  1)  x  (Mx  +  1)  can  be  formed  such  that 


G=[g0        gl         •••       gMx 


(4.12) 
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Now  the  second  step  of  the  row-column  decomposition  can  be  performed.    From 
Eq.   (4.5),  for  a  particular  wlfc, 

U  x 


0(wifc,n2)=  ^T  a(n1,n2)coso;lfcn1. 


(4.13) 


n,  =0 


For  a  particular  n2 ,  a  system  of  linear  equations  can  be  formed  by  evaluating  Eq. 
(4.13)  for  each  k,  k  =  0, 1, . . . ,  Mi .  In  matrix  form,  this  becomes 


g{vio,n2) 

?(wn,n2) 


1       cosWio 

1  COSCJn 

1     cosw1Ml 


cos  Mi  uji  1 
cosM1a;1Ml  . 


a(0,n2) 
a(l,n2) 


(4.14) 


Equation  (4.14)  can  be  expressed  in  vector-matrix  notation  as 


Ba, 


(4.15) 


where  rn3  is  a  (Mx  +  1)  x  1  column  vector  formed  from  the  elements  of  the  n2th 
row  of  G,  B  is  a  (Mi  +  1)  x  (Mx  +  l)  square  matrix  and  a  is  a  (Mx  +  1)  x  1  column 
vector.  Thus  a  is  found  from  the  relation 


a  =  B1T 


n-i  • 


(4.16) 


The  vector  a  represents  the  values  of  a(rii ,  n2)  for  a  specified  value  of  n2 .  To  obtain 
the  entire  (Mi  +  1)  x  (M2  +  l)  sequence  a(ni,n2),  the  process  described  by  Eq. 
(4.16)  is  repeated  for  each  value  of  n2,  n2  =  0,1,..., M2.  Therefore,  a  total  of 
(M2  +  1)  systems  containing  (Mx  +  1)  linear  equations  each  are  involved  in  the 
second  step  of  the  method.  Finally,  the  impulse  response  of  the  resulting  filter  can 
be  evaluated  using  the  relationship  of  Eq.  (4.2). 

The  following  design  example  illustrates  how  this  method  can  be  employed. 
For  this  example,  a  2-D  circularly  symmetric  lowpass  filter  characteristic  is  desired 
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with  unity  gain  in  the  passband  region  {(ct>i,k;2)  :.0  <  \fu\  +  w2  <  0.27r}  and  zero 


gain  in  the  stopband  region  {(w1,w2)  :  0.47T  <  ywj  +  ^2  <  ?r}-  The  region  of 
support  for  the  impulse  response  is  specified  to  be  17  x  17.  The  arrangement  of 
frequency  samples  shown  in  Figure  4.2  is  to  be  used.  Note  that  the  sampling  is 
not  uniform.  When  samples  are  placed  in  the  transition  region,  a  linear  slope 
characteristic  between  the  two  bands  is  used  in  specifying  the  desired  frequency 
response  value  at  those  samples,  as  was  done  in  the  uniform  sampling  approach. 
Contour  and  perspective  plots  of  the  frequency  response  for  the  resulting  filter  are 
shown  in  Figure  4.3. 


2.1 


wl 


Figure  4.2.  Arrangement  of  Samples  for  17  x  17  Lowpass  Filter. 

One  rather  obvious  extension  of  the  nonuniform  sampling  method  involves 
again  specifying  an  Ny  x  N2  impulse  response  region  of  support,  but  in  this  case 
performing  the  interpolation  in  the  first  step  of  the  process  on  wt .  Suppose  any 
(M2  +  1)  distinct  values  of  u>2  are  selected  on  the  interval  0  <  ui2  <  tt-  These  are 
denoted  as  {u2i,  I  —  0, 1, . . . ,  A/2}-    For  each  value  of  u2l,  any  [Mx  +  l)  distinct 
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(a) 


1.0  } 


00 
Figure  4.3.  Frequency  Response  of  17  x  17  Lowpass  Filter. 
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values  of  Wj  are  chosen  on  the  interval  0  <  u^  <n.  These  values  are  denoted  as 
{ulk(l),  k  =  0, 1, . . . ,  Mi  } .  Such  an  arrangement  of  selected  samples  is  illustrated 
in  Figure  4.4  for  the  case  A/x  =  A/2  =  2.  Now  suppose  H(u)l,u2)  is  specified  at 
each  selected  sample  (wlfc  (/),w2,)  for  k  —  0, 1, . . . ,  Ai\  and  /  =  0, 1,. . .  ,  M2. 


CJ, 


—  -♦ 


0 


•  - 


>io(l)        wu(l)       w12(l) 


^io(O) 


<-:i(0)     u,12(0) 


7T       OJ. 


Figure  4.4.  Alternate  Nonuniform  Sampling  Approach. 
Then  there  is  a  unique  solution  for  a[nl  ,  n2)  which  can  be  obtained  from 


H[u}lk(l),u}2i)  =  2J    zZ  aini,n2)cosulk(l)n1cosuj2tn2; 


0  <  /c  <  Mx ,     0  <  /  <  M2 . 


(4.17) 


The  approach  to  the  solution  of  the  above  problem  and  proof  of  its  existence  and 
uniqueness  can  be  performed  in  a  straightforward  manner  analogous  to  that  of  the 
basic  approach. 
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It  is  worthy  to  note  that  the  sampling  geometry  required  by  the  uniform 
sampling  method  of  Chapter  II  is,  in  reality,  a  special  case  of  the  general  nonuni- 
form sampling  form  discussed  here.  Therefore,  the  results  of  the  uniform  sampling 
method  can  always  be  achieved  by  constraining  the  location  of  samples  to  the  uni- 
form Cartesian  grid.  However,  the  method  of  solution  is  less  efficient  than  the 
IDFT  approach  and  so  would  probably  not  be  used  for  this  particular  case. 

B.     EXISTENCE  AND  UNIQUENESS 

A  significant  aspect  of  the  nonuniform  sampling  method  is  that  a  unique  solu- 
tion matching  the  specified  frequency  response  at  each  sample  point  is  guaranteed, 
provided  samples  are  located  as  described.  Proof  of  existence  and  uniqueness  of  a 
sequence  h(n1,n2)  which  satisfies  the  interpolating  conditions  follows. 

The  first  step  in  the  nonuniform  sampling  approach  involves  solution  of  [My  + 
1)  systems  of  linear  equations  of  the  form  of  Eq.  (4.9)  or  Eq.  (4.10).  A  unique 
solution  exists  for  each  of  these  systems.  Recall  that  Eq.  (4.9)  resulted  from 
applying  the  interpolating  conditions  to  Eq.  (4.6).  Through  use  of  the  Chebyshev 
polynomials,  Eq.  (4.6)  can  be  expressed  as 

Mi 

#(o>i,w3)=   J^  c(w1,n2)(cosw2)ns,  (4.18) 

nj  =0 

where  a  unique  relation  exists  between  the  coefficients  g(u>l,n2)  and  c(oj1,n2). 
Through  the  subtitution  x  =  cos  u>2 ,  the  above  becomes 

Mi 

H(ux,x)  =   ^2  c(w!,n2)xn3.  (4.19) 

n7-Q 

The  set  {l,x,x2, . . .  ,xM3}  is  known  to  be  unisolvent  on  any  interval  [a,b]  [Ref.  9: 
p.  31].  Therefore,  a  unique  interpolation  polynomial  of  the  form  of  Eq.  (4.19)  can 
be  found  given  any  (M2  +  1)  distinct  values  of  x.  From  the  relation  xt  =  cos^;, 
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as  long  as  distinct  w2l{k)  samples  are  specified  on  the  interval  {0  <  u>2  <  7r},  the 
corresponding  xt  will  be  distinct.  Due  to  the  unique  relationship  between  Eqs. 
(4.6),  (4.18),  and  (4.19),  the  interpolation  implied  by  Eq.  (4.10)  will  yield  a  unique 
solution  for  gfc,  k  =  0,1,...,  Mx.  This  implies  that  the  matrix  G  defined  in  Eq. 
(4.12)  will  be  unique. 

The  second  step  of  the  nonuniform  sampling  approach  involves  solution  of 
(M2  +  1)  systems  of  linear  equations  of  the  form  of  Eq.  (4.15).  Each  row  of  matrix 
B  in  Eq.  (4.15)  forms  a  vector  of  Chebyshev  polynomials  which  form  a  unisolvent 
set.  Therefore,  B  is  invertible.  Since  a  unique  matrix  G  was  found  in  the  first  step, 
the  column  vectors  Tn2  are  unique  for  each  n2  —  0, 1, . . . ,  M2 .  Therefore,  a  unique 
solution  for  a^  ,n2)  exists  which  satisfies  the  overall  interpolation  problem.  From 
the  relationship  indicated  by  Eq.  (4.2),  this  implies  existence  and  uniqueness  for 
the  resulting  impulse  response  h(n1,n2). 

Ill-conditioned  results  normally  occur  in  interpolation  problems  involving  one 
variable  when  one  sample  is  located  in  close  proximity  to  another.  Since  the 
nonuniform  sampling  method  uses  a  two-step  process  involving  one-dimensional 
interpolating  functions,  the  same  result  can  be  expected  to  hold  here.  Provided 
samples  are  not  placed  too  close  together,  the  resulting  filter  design  should  not  be 
ill-conditioned. 

C.     COMPUTATIONS 

While  the  nonuniform  sampling  method  retains  less  flexibility  than  the  ar- 
bitrary sampling  method  of  Chapter  III,  it  is  more  computationally  efficient  for 
filters  of  the  same  order.  As  in  the  analysis  for  the  arbitrary  sampling  case,  in- 
verting an  N  x  N  matrix  is  assumed  to  require  on  the  order  of  iV3  floating  point 
operations  and  additions,  and  inversion  of  the  matrix  is  assumed  to  dominate  the 
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computations  required  in  solving  a  system  of  linear  equations.  A  brief  analysis  of 
the  number  of  operations  required  for  the  nonuniform  sampling  method  follows. 

For  a  filter  with  a  specified  impulse  response  region  of  support  of  size  Ni  xJV2, 
step  one  of  the  design  process  requires  solution  of  a  total  of  (Mx  +  1)  systems  of 
(M2  +  1)  linear  equations  each.  The  number  of  floating  point  operations  involved 
in  this  step  is  therefore  on  the  order  of  (Mt  +  l)(Af2  +  l)3. 

Step  two  of  the  design  process  requires  solution  of  a  total  of  (M2  +  1)  systems 
of  (Mi  + 1)  linear  equations  each.  The  number  of  required  floating  point  operations 
in  this  step,  then,  is  on  the  order  of  (M2  +  l)(-Mi  +  l)3. 

The  total  number  of  floating  point  operations  for  the  entire  nonuniform  frequency- 
sampling  2-D  FIR  filter  design  method  is  the  combined  amount  from  the  above 
two  steps.  The  total  number,  therefore,  is  on  the  order  of 

Cnu  a  (Mx  +  1)(M2  +  l)3  +  (Mx  +  1)3(M2  +  1) 

JVi  -  1        \/N2-l        \3       (Nx-\        \ 3  (N2  -  1 

+  l)f-~ —  +  1)    +     — ~z — +1]    l-^-z —  +  1 


1  "  (4-20) 

=  jjW  +  l)(N2  +  l)3  +  -(Nx  +  1Y(N2  +  1) 

For  the  particular  case  in  which  Nx  =  N2  =  N,  Eq.  (4.20)  reduces  to 

CnucijV4.  (4.21) 

o 

Thus,  the  nonuniform  sampling  method  is  much  more  computationally  effi- 
cient than  the  arbitrary  sampling  method,  yet  is  still  quite  flexible  while  providing 
guarantee  of  a  unique  design  solution. 
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V.     A  COMPARISON  OF  FREQUENCY-SAMPLING  DESIGN  METHODS 

THROUGH  DESIGN  EXAMPLES 

Up  to  this  point  the  theory  behind  three  fundamentally  different  approaches 
to  the  frequency-sampling  design  of  2-D  FIR  filters  has  been  developed.  With 
knowledge  of  these  approaches,  it  is  of  interest  to  investigate  the  results  of  applying 
each  of  the  methods  to  some  practical  design  problems.  In  this  chapter,  some  design 
problems  are  formulated  and  the  results  from  each  of  the  methods  are  compared 
in  terms  of  performance,  computations,  and  flexibility. 

A.      APPROACH 

In  order  to  compare  the  three  basic  design  methods  discussed  in  this  thesis, 
three  different  design  problems  are  described.  To  achieve  a  valid  comparison,  sev- 
eral parameters  were  fixed  for  each  design.  The  region  of  support  for  the  filter 
impulse  response  in  all  cases  was  specified  to  be  17  x  17.  A  square  region  of  sup- 
port is  most  commonly  used  in  2-D  filter  design  and  avoids  some  of  the  problems 
discussed  in  the  latter  part  of  Chapter  IV.  The  specified  dimensionality  of  the  im- 
pulse response  is  large  enough  to  ensure  reasonable  results  but  not  so  large  as  to 
make  implementation  impractical. 

Each  design  method  requires  specification  of  a  desired  frequency  response 
Hd{ojx  ,oj2)-  For  each  design  method,  ifd(wi  ,w2)  is  fixed,  with  unity  gain  through- 
out the  passband,  zero  gain  throughout  the  stopband,  and  a  linear  slope  charac- 
teristic in  the  transition  band(s).  The  frequency  response  value  assigned  to  each 
sample  point  in  the  two-dimensional  frequency  domain  is  the  value  of  Hd(u1 , w2) 
at  that  point. 
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For  each  design  problem,  circularly  symmetric  filters  are  prescribed.  Although 
more  general  symmetries  are  allowable  for  each  method,  circularly  symmetric  de- 
signs adequately  illustrate  the  salient  characteristics  of  the  methods. 

A  few  other  points  are  of  note.  For  these  design  examples,  optimimization 
techniques  were  not  utilized  in  locating  frequency  samples  for  any  of  the  methods. 
Frequency  sample  locations  were  selected  strictly  on  a  trial-and-error  basis  by  the 
human  designer  for  the  arbitrary  and  nonuniform  sampling  methods.  The  resulting 
frequency  responses  of  the  filters  are  compared  qualitatively  in  terms  of  passband 
and  stopband  ripple,  transition  band  width,  and  circularity.  No  quantitative  com- 
parisons are  made  in  terms  of  best  equi-ripple  approximation  to  an  ideal  response 
characteristic,  etc.  It  is  well  known  that  a  better  quality  filter  can  usually  be  ob- 
tained by  increasing  the  number  of  terms  in  the  impulse  response  (i.e.  increasing 
the  filter  order).  In  comparing  frequency-sampling  methods  with  the  same  support 
for  the  impulse  response  in  each  case,  substantially  better  filters  are  not  expected 
to  result  from  one  particular  method;  rather,  it  is  of  interest  to  investigate  how 
design  method  flexibility  can  be  used  to  enhance  certain  desirable  characteristics 
in  the  resulting  filter  frequency  response. 

B.      DESIGN  EXAMPLES 

Three  different  design  problems  are  now  presented.  The  three  frequency- 
sampling  design  methods  are  applied  to  each  and  a  comparative  analysis  of  the 
resulting  filter  frequency  responses  is  performed. 

1.     Lowpass  Filter  No.  1 

For  this  example,  a  circularly  symmetric  2-D  lowpass  filter  frequency  re- 


sponse is  desired  such  that  it  has  unity  gain  in  the  passband  region  {0  <  \fw1  +  oj%  < 


0.47r}  and  zero  gain  in  the  stopband  region  {0.67T  <  \J(jj\  +  u>%  <  7r}.  First,  a  17x  17 
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filter  was  designed  to  match  the  above  characteristics  using  the  uniform  sampling 
method.  Contour  and  perspective  plots  of  the  filter  frequency  response  are  shown 
in  Figure  5.1.  Next,  a  17  x  17  filter  was  designed  using  the  nonuniform  sampling 
method  with  the  arrangement  of  samples  shown  in  Figure  5.2.  The  resulting  fil- 
ter frequency  response  is  illustrated  in  Figure  5.3.  Finally,  the  arbitrary  sampling 
method  was  applied  to  this  design  problem.  The  selection  of  samples  was  as  shown 
in  Figure  5.4,  and  the  frequency  response  of  the  resulting  filter  can  be  seen  in 
Figure  5.5. 

In  comparing  the  results  of  the  three  methods,  the  quality  of  the  result- 
ing filters  varies  slightly  in  each  case  but  the  arbitrary  and  nonuniform  sampling 
designs  appear  no  better  than  that  of  the  uniform  sampling  case.  The  uniform 
sampling  method  required  only  2  x  173  =  9826  floating  point  operations,  while 
the  nonuniform  sampling  method  required  on  the  order  of  2  x  94  ~  13,000  float- 
ing point  operations.  The  arbitrary  sampling  method  required  on  the  order  of 
1/64  x  176  ~  380,000  floating  point  operations.  In  looking  back  at  the  desired 
characteristics  of  this  filter,  the  passsband  and  stopband  regions  comprise  fairly 
equal  areas  of  the  (wi,^)  plane.  Sampling  on  a  uniform  Cartesian  grid  allows 
an  adequate  number  of  samples  in  all  bands  and  hence  little  seems  to  be  gained 
through  nonuniform  spacing  of  samples. 
2.     Lowpass  Filter  No.  2 

For  the  second  example,  another  circularly  symmetric  2-D  lowpass  filter 
characteristic  is  desired.  In  this  case,  the  filter  is  to  have  unity  gain  in  the  passband 


{0  <  \/w|+q|  <  0.27r}  and  zero  gain  in  the  stopband  {0.47T  <  \Ju>\  +  u>l  <  7r}. 
This  particular  example  was  selected  since  it  illustrates  a  case  in  which  the  desired 
passband  encompasses  a  much  smaller  portion  of  the  region  K  of  the  (u^  ,  u>2)  plane 
than  the  stopband. 
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Figure  5.1.  Frequency  Response  of  Lowpass  Filter  Design  No.  1, 
Uniform  Sampling  Method. 
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Figure  5.2.  Arrangement  of  Samples  for  Lowpass  Filter 
Design  No.  1,  Nonuniform  Sampling  Method. 
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Figure  5.3.  Frequency  Response  of  Lowpass  Filter  Design  No.  1, 
Nonuniform  Sampling  Method. 


53 


•  •  • 


in 


a  ^ 


• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

- 

• 

• 

• 

• 

• 

• 

• 

- 

• 

• 

• 

• 

■ 

_ 

• 

♦ 

• 

■ 

• 

m 

• 

• 

■ 

• 

• 

• 

• 

• 

• 

• 

1  ♦ 

1 

•    1 

1 

. 

\ 

• 

• 

« 

• 

0.5 


1.0 


1.5 
u1 


2.0 


2.5 


3.0 


Figure  5.4.  Arrangement  of  Samples  for  Lowpass  Filter 
Design  No.  1,  Arbitrary  Sampling  Method. 
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Figure  5.5.  Frequency  Response  of  Lowpass  Filter  Design  No.  1, 
Arbitrary  Sampling  Method. 


55 


The  use  of  the  uniform  sampling  method  allows  the  placement  of  only  four 
samples  in  the  passband  region.  The  resulting  filter  frequency  response,  shown  in 
Figure  5.6,  has  a  rather  smooth  passband  and  stopband  characteristic,  but  the 
passband  is  smaller  in  area  than  desired  and  can  be  seen  to  have  a  rather  squared 
shape  in  the  contour  plot. 

The  specified  selection  of  frequency  samples  for  the  nonuniform  sampling 
method  is  illustrated  in  Figure  5.7.  Note  that  five  samples  are  now  placed  within 
the  desired  passband  region.  Figure  5.8  represents  the  frequency  response  of  the 
resulting  filter.  In  this  case,  the  size  of  the  resulting  passband  more  closely  matches 
the  desired  characteristic  and  the  passband  is  clearly  more  circular  in  shape.  These 
characteristics  are  gained,  however,  at  the  expense  of  introducing  a  greater  amount 
of  ripple  in  the  stopband. 

Figure  5.9  shows  a  chosen  arrangement  of  samples  for  the  arbitrary  sam- 
pling method.  In  this  instance,  six  samples  are  now  contained  in  the  specified 
passband  region  while  no  samples  are  placed  in  the  transition  region.  The  result- 
ing filter  frequency  response  is  illustrated  in  Figure  5.10. 

Note  that  for  the  latter  result,  the  passband  is  fairly  flat  and  quite  circular, 
and  the  roll-off  is  quite  sharp.  This  occurs  since  the  flexibility  of  the  design  allows 
frequency  samples  to  be  placed  along  the  perimeters  of  the  specified  passband 
and  stopband  regions  for  the  latter  two  design  methods.  Because  of  this,  the  two 
methods  which  allow  flexibility  in  sampling  can  be  used  to  provide  greater  control 
over  the  resultant  width  of  the  transition  band.  Of  course,  this  is  achieved  at  the 
expense  of  increased  ripple.  Since  the  arbitrary  and  nonuniform  sampling  methods 
do  not  require  sampling  on  a  Cartesian  grid,  they  are  generally  better  able  to 
conform  to  non-rectangularly  shaped  contours  such  as  circles. 
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Figure  5.6.  Frequency  Response  of  Lowpass  Filter  Design  No.  2, 
Uniform  Sampling  Method. 
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Figure  5.7.  Arrangement  of  Samples  for  Lowpass  Filter 
Design  No.  2,  Nonuniform  Sampling  Method. 
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Figure  5.8.  Frequency  Response  of  Lowpass  Filter  Design  No.  2, 
Nonuniform  Sampling  Method. 
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Figure  5.9.  Arrangement  of  Samples  for  Lowpass  Filter 
Design  No.  2,  Arbitrary  Sampling  Method. 
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Figure  5.10.  Frequency  Response  of  Lowpass  Filter  Design  No.  2, 
Arbitrary  Sampling  Method. 
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3.     Bandpass  Filter 

The  final  design  example  illustrates  how  each  of  the  frequency-sampling 
techniques  can  be  used  in  the  design  of  a  2-D  circularly  symmetric  bandpass  filter. 


Here,  the  desired  filter  has  unity  gain  in  the  passband  region  {0.47T  <  yjuj\  +  w\  < 


0.67t},  and  zero  gain  in  the  stopbands  defined  by  {0  <  yjw\  +  u>\  <  0.2iv}  and 
{0.87T  <  \Joj1  +  u>2  <  tt}-  The  filter  resulting  from  application  of  the  uniform 
sampling  method  has  the  frequency  response  illustrated  in  Figure  5.11.  For  the 
nonuniform  sampling  method,  the  arrangement  of  samples  shown  in  Figure  5.12 
was  used.  The  resulting  filter  frequency  response  is  shown  in  Figure  5.13.  Figure 
5.14  illustrates  the  arrangement  of  samples  selected  for  the  arbitrary  sampling 
method.  The  resulting  filter  has  the  frequency  response  of  Figure  5.15. 

In  comparing  these  results,  the  nonuniform  sampling  method  produces  a 
somewhat  flatter  and  wider  passband  than  the  uniform  sampling  method,  but  re- 
sults in  some  notable  ripples  in  the  outer  stopband  which  are  not  found  in  the 
uniform  sampling  case.  The  passband  for  the  arbitrary  sampling  design  is  even 
more  wide  and  flat,  but  at  least  for  this  particular  selection  of  samples,  has  unac- 
ceptably  high  ripple  in  both  stopband  regions. 

The  point  to  be  made  in  comparing  these  methods  is  not  that  the  nonuni- 
form or  arbitrary  sampling  methods  necessarily  produce  better  filters,  but  rather 
that  they  possess  more  flexibility  and  hence  can  be  used  to  enhance  a  particular 
desirable  feature  at  the  expense  of  another.  The  advantage  gained  through  design 
flexibility  may  well  outweigh  the  cost  of  additional  computations  for  many  design 
problems. 
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Figure  5.11.  Frequency  Response  of  Bandpass  Filter  Design, 
Uniform  Sampling  Method. 
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Figure  5.12.  Arrangement  of  Samples  for  Bandpass  Filter 
Design,  Nonuniform  Sampling  Method. 
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Figure  5.13.  Frequency  Response  of  Bandpass  Filter  Design, 
Nonuniform  Sampling  Method. 
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Figure  5.14.  Arrangement  of  Samples  for  Bandpass  Filter 
Design,  Arbitrary  Sampling  Method. 
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Figure  5.15.  Frequency  Response  of  Bandpass  Filter  Design, 
Arbitrary  Sampling  Method. 
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VI.     CONCLUSIONS 

In  this  thesis,  various  approaches  to  the  frequency-sampling  design  of  two- 
dimensional  FIR  digital  filters  were  investigated.  The  traditional  frequency-sampling 
approach,  referred  to  as  the  uniform  sampling  method  here,  involves  specifying  the 
frequency  response  at  sample  points  located  at  the  vertices  of  a  uniform  Cartesian 
grid  in  the  (oj1,u2)  plane.  Through  use  of  the  IDFT,  the  method  is  quite  compu- 
tationally efficient,  and  can  be  made  even  more  efficient  if  an  FFT  is  used.  The 
method  possesses  no  flexibility  in  matching  design  parameters,  however.  Another 
frequency-sampling  approach  arises  from  specifying  the  frequency  response  at  ar- 
bitrarily selected  distinct  samples  in  the  (u^,^)  plane  and  then  solving  a  large 
system  of  linear  equations.  This  method,  referred  to  as  the  arbitrary  sampling 
method  in  this  thesis,  possesses  great  flexibility  but  is  computationally  intensive 
and  prone  to  the  occurrence  of  degeneracies  in  which  a  unique  solution  does  not  ex- 
ist. Finally,  a  new  frequency-sampling  approach,  termed  the  nonuniform  sampling 
method,  was  introduced.  This  method  somewhat  constrains  the  location  of  fre- 
quency samples,  but  still  provides  some  free  parameters  and  is  reasonably  efficient 
in  computation. 

The  results  of  this  thesis  reflect  how  computational  efficiency  and  design  flexi- 
bilty  relate  to  the  number  of  constraints  on  sample  location.  The  methods  described 
here  which  impose  more  constraints  are  more  computationally  efficient.  Methods 
with  less  constraints  were  found  to  be  more  flexible  in  matching  a  desired  frequency 
response. 

A  problem  posed  by  both  the  arbitrary  and  nonuniform  sampling  methods  is 
how  to  properly  select  the  samples  in  order  to  achieve  an  acceptable  filter.  Various 
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algorithms  for  the  design  of  2-D  optimal  filters  have  been  developed  which  involve 
the  solution  of  linear  equations,  analogous  to  those  for  the  arbitrary  sampling 
method,  in  an  iterative  approximation  process.  Although  they  have  been  shown  to 
converge,  the  algorithms  necessarily  take  special  measures  to  deal  with  degeneracies 
and  are  very  computationally  intensive.  [Refs.  9,  10] 

The  nonuniform  sampling  method  warrants  a  closer  look  because  of  its  rela- 
tive computational  efficiency.  Due  to  constraints  on  the  sampling  geometry,  the 
nonuniform  sampling  method  is  probably  not  well-suited  to  the  design  of  globally 
optimal  filters.  If  an  efficient  algorithm  can  be  developed  which  will  locate  a  set  of 
points  satisfying  the  non-uniform  method  constraints,  and  if  this  set  guarantees  a 
filter  which,  while  not  necessarily  optimal,  nevertheless  is  of  acceptable  quality,  the 
non-uniform  sampling  method  developed  in  this  thesis  may  prove  to  be  of  great 
significance. 
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APPENDIX 

A  2-D  FIR  FILTER  DESIGN  BASED  ON  A  BIVARIATE  EXTENSION 
OF  NEWTON'S  POLYNOMIAL  INTERPOLATION  FORMULA 

In  the  course  of  this  thesis  research,  an  extension  of  Newton's  polynomial 
interpolation  method  in  two  variables  was  studied  and  applied  to  the  design  of  two- 
dimensional  digital  filters.  Such  an  extension  of  Newton's  method  was  proposed  by 
Diamessis  [Ref.  11].  Its  application  as  a  2-D  FIR  filter  frequency-sampling  design 
method  was  presented  at  the  ICASSP  87  conference  [Ref.  12]  and  is  summarized 
here. 

A.     2-D  NEWTON  INTERPOLATION 

The  bivariate  Newton  interpolation  problem  for  the  2-D  quadratic  case  can 
be  described  as  follows.  Given 

1)  a  set  of  points  in  the  (x,y)  plane: 

S/  ={{xo,yo),(x0,y1)i(x0,y2),(x1,y0),(x1,y1),(x2,y0)}  (A  -  1) 

2)  a  set  of  values  at  those  points: 

t>   A     I    /(X0i  Vo)   =  /oo,        f[Xo,yi)   =  fo\  >        f{Xo,y2)   =   /o2,    I  ( A     -  9^ 

"l/(*i,l(b)4/iQ,     /(*i,Vi)4/ii,     /(x2,y0)^/20  J  [  ] 

3)  a  quadratic  2-D  polynomial: 

p{x,y)  =  a00  +  al0x  +  a01y  +  a20x2  +  axlxy  +  a02y2  (A  -  3) 

find 

a  —  [a0o  i  flio  )aoii  a2o  i  an  5  ao2  ]  (A  —  4) 
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such  that  the  following  interpolating  conditions  are  satisfied 


p(zo>!/o)  =  /oo,     p(so,yi)  = /on     P(zo,!/2)  =  /02 

p(zi,y0)  =  /m,    p(xi,yj  =  fn,    p{x3iy0)  =  /20 


If  the  polynomial  of  Eq.  (A-3)  is  rewritten  in  the  form 
p(x,y)  =  c00  -f  c10(i  -  x0)  +  c01(y  -  y()) 

+  c20  (x  -  Xo  )  (x  -  xx )  +  ci ,  (x  -  Xu  )  (y  -  y„  ) 

+  c02(y  -  y0)(y  -  yx), 
then  the  interpolation  problem  involves  finding  the  coefficient  vector 


(A -5) 


(A -6) 


A     f  ]T 

C    —     Coo,Cio,Coi,C2o,C11,C(j2 


(A -7) 


instead  of  a.   The  system  of  equations  formed  by  application  of  Eq.    (A-2)  to  Eq. 
(A-6)  is 


1            0 

0 

1        {x1   -  X0) 

0 

1             0 

[yi  -  yo) 

1      (i2  -  I0) 

0 

1      (ii  -  x0) 

(i/i  -  yo) 

1             0 

(y2  -  yo) 

0  0 

0  0 

0  0 

(x2-  z0)(x2  -  Xi)  0 

0  (xi  -  xo)(yi  ~  yo] 
0  0 


"/oo 
/io 

/oi 
/20 

Ai 

/02 


(y2  -  yo)(y2  -  yi). 


coo 

c10 
C()l 

C20 

Cll 

Lc02 


(A -8) 


or,  in  vector-matrix  notation, 


NonC    =f. 


(A -9) 
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Since  N2D  is  lower  triangular,  c  can  be  solved  for  recursively.  Also,  N2Z)  is 
non-singular  provided  the  x{  are  distinct  and  the  t/;  are  distinct.  Hence  a  unique 
solution  for  c  can  always  be  found.  The  coefficient  vector  a  is  related  to  c  by  the 
relation 

c  =  U2Da,  (A -10) 

where  U2D  is  an  upper  triangular  matrix  which  transforms  the  Newton  coefficient 
vector  c  into  the  coefficient  vector  a.  Therefore,  a  can  be  found  through  solution 
of  the  two  triangular  systems  indicated  by  Eqs.  (A-9)  and  (A-10). 

B.     NEWTON  FILTER  DESIGN  METHOD 

The  selection  of  frequency  samples  for  the  Newton  design  method  follows  from 
the  two  variable  Newton  interpolation  method  discussed  above.  This  approach  is 
a  2-D  extension  of  a  one-dimensional  filter  design  proposed  by  Schuessler  [Ref.  13: 
p.  257].  If  a  square  region  of  support  on  the  first  quadrant  of  extent  N  x  N  is 
initially  specified  for  the  filter  impulse  response,  the  frequency  response  of  such  a 
2-D  FIR  filter  with  linear  phase  can  be  expressed  as 

N  -1    N  -  1 

H(ultua)=J2    Y,  h{n1,n2)e-i»^e->">n>.  (A  -  11) 

A  total  of  (M  +  1)(M  +  2)/2  frequency  samples  are  then  chosen  on  region  K  of  the 
[u>! ,  w2)  plane,  where  M  —  (N  —  l)/2  and  K  —  {{w\ ,  w2),0  <  uj1  <  7r;0  <  u>2  <  7r}. 
The  selected  samples  are  of  the  form 

{((wi^.WaO.Afe  =0,1,. ...A*),*!  =0,1,.  ..,N}.  (A  -  12) 

A  support  set  of  the  above  form  is  referred  to  as  triangular.  Note  that  the  samples 
need  not  be  uniformly  spaced.  Since  there  is  no  constraint  such  that  w,-0  <  u>il  < 
u>{j  •  •  •,  the  sampling  geometry  is  not  necessarily  of  a  triangular  shape  on  region 
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K.  Two  such  support  sets  are  shown  in  Figure  A.l.  Given  a  specified  value  for 
H(u>ik  ,u)2k)  at  each  sample  location,  the  above  support  set  ensures  existence  of  a 
unique  solution  to  the  design  problem. 

The  Newton  design  imposes  the  quadrant  symmetry  condition 

h(ni,n2)  =  h(N  —  1  —  nlin2)  =  /i(n.1}iV  —  1  —  n2).  (A  —  13) 

Application  of  Eq.  (A-13)  to  Eq.  (A-ll)  results  in  the  expression 

M  -  1 

H0{u1,u2)=h{M,M)  +  2  Y^  h(M,n2)TM  _n,  (cos u2) 

(A -14) 

M  -  1     2M  v  ' 

+  2  22    5Z  h,{ni^n2)TM-ni(cosu1)TM-n3[cosu2)i 

where  Tn  (x)  is  the  nth  degree  Chebyshev  polynomial  in  x  and  H0(oj1,tL>2)  is  the 
zero-phase  frequency  response  related  to  H{wl,w2)  by 

H{ul9u2)  =  e-Ju*n*e-'u'n'H0{<jjltu2).  (A  -  15) 

Now,  if  the  impulse  response  region  of  support  is  restricted  such  that 

M»i,n2)=0,      lni  ~  M\  +  ln2  -  M\  >  M,  (A  -  16) 

and  if  the  transformations 

cosu^  =  1  —  2x     and      cosu^  =  1  —  2y  (A  —  17) 

are  applied,  then  the  zero-phase  frequency  response  of  Eq.  (A- 14)  is  mapped  into 
a  bivariate  polynomial  of  the  Newton  form 

M  t -  1  M  j-1 

H{x,y)  =c00  +^c,0  JJ(x-xp)  +  J^  coy  JI(y-yg) 
1-1         P  =  o  y.,         ,=  o 

M     i : -  1  i-j-1  j-1 

+^2^2c<-jj  n  (x_xp) n>~^)- 

t=2  J  =  l  p=0  q=0 
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Figure  A.l.  Possible  Sample  Arrangements  for  Bivariate  Newton 
Interpolation,  Quadratic  Case. 
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The  design  problem  involves  first  mapping  the  specified  frequency  samples 

(u>lfc ,  u;2fc )  into  samples  (xk  ,  yk )  through  the  relations  of  Eq.  (A-17),  and  then  using 

the  Newton  interpolation  method  previously  described  to  find  the  coefficients  ci2 . 

By  applying  Eq.  (A-17)  to  Eq.  (A-18),  and  through  the  substitution 

e?u*  +e-JW-         Zi  +  z~l 
2  2 

a  filter  structure  of  the  form 


H(Zl,z2)=c00z;"z;"  +^c,o2l-(M-"22-w  JJ-i(i-Mlp<-i+^-») 

t=l  p=0 

M  j  -  1 

+x:co^rAf^,"-j)n-i(i-^^i+^2) 


4 

3=1  <J  =  0 


M     i-  1 

f  -(M-i+j) 

'2 


t  =  2  j  =  1 

» -  j  - 1     -  y  - 1 


p=0  g=0 

(A  -  19) 
results,  where  yulp  =  2cosu;1p  and  fi2q  —  2coso;2q- 

C.     DISCUSSION 

The  advantages  hoped  to  be  gained  by  design  of  the  Newton  filter  were  com- 
putational efficiency  through  the  recursive  nature  of  the  solution  and  use  of  the 
permanence  property  [Ref.  11].  The  latter  property  allows  the  specification  of 
additional  sample  values  for  an  interpolation  problem  without  having  to  redeter- 
mine previously  computed  coefficients.  These  advantages,  however,  apply  only  if 
the  filter  is  implemented  using  the  structure  implied  by  Eq.  (A-19),  where  the  co- 
efficients c{j  represent  the  gains.  Such  a  structure,  however,  requires  more  delays 
than  a  direct  form  structure  of  the  same  order  which  uses  the  impulse  response 
coefficients  h(nx,n2)  for  the  filter  gains,  and  hence  is  not  very  attractive. 
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An  alternate  procedure  might  be  to  use  the  Newton  method  to  find  the  im- 
pulse response  coefficients  and  implement  the  filter  in  direct  form.  Obtaining  the 
impulse  response  coefficients  of  the  filter  from  the  set  of  coefficients  c,y  requires 
solution  of  two  triangular  system  analogous  to  those  of  Eqs.  (A-9)  and  (A-10). 
Carrying  out  this  two-step  process,  solving  first  for  c  and  then  for  h(nl,n2)  to 
get  a  direct  form  implementation,  is  equivalent  to  solving  a  system  of  equations 
via  LU-decomposition.  This  particular  system  is  the  system  of  linear  equations 
resulting  from  the  direct  application  of  the  interpolating  conditions  to  Eq.  (A-ll), 
where  the  limits  of  the  summation  in  Eq.  (A-ll)  are  altered  to  reflect  the  Newton 
filter  region  of  support.  This  is  the  same  approach  used  in  the  arbitrary  sam- 
pling method  of  Chapter  III,  if  the  matrix  inversion  is  performed  through  use  of 
LU-decomposition,  and  therefore  cannot  be  considered  any  more  computationally 
efficient. 

The  Newton  method  may  be  less  significant  in  terms  of  computational  effi- 
ciency or  practical  filter  implementation  than  it  is  for  studying  the  formulation 
of  the  two-dimensional  frequency-sampling  design  problem.  In  Chapter  III,  it  was 
shown  that  existence  and  uniqueness  of  a  solution  to  the  interpolation  problem 
for  arbitrarily  placed  samples  and  a  fixed  impulse  response  region  of  support  was 
not  guaranteed  for  the  two-dimensional  case.  With  specification  of  a  triangular 
support  set  and  a  polynomial  of  the  form  of  Eq.  (A-6),  the  Newton  method  pro- 
vided a  unique  solution.  If  another,  less  constrained  choice  of  samples  is  specified 
and  if  the  form  of  the  polynomial  is  altered  to  include  basis  functions  such  that 
the  interpolating  problem  becomes  a  lower  triangular  system  of  equations  akin  to 
Eq.  (A-8)  with  a  non-zero  main  diagonal,  then  a  unique  interpolating  polynomial 
of  the  new  form  will  exist.  Since  the  specified  basis  functions  for  the  polynomial 
interpolation  problem  directly  relate  to  the  region  of  support  for  the  impulse  re- 
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sponse  in  the  filter  design  problem,  the  Newton  method  can  be  used  to  illustrate 
a  connection  between  the  particular  sampling  geometry  and  an  impulse  response 
region  of  support  which  will  ensure  a  unique  solution. 
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